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Abstract 


I  propose  to  analyze  incompressible  flow  of  a  Newtonian 
Fluid  past  a  vertical,  flat  plate  with  thermal  and  magnetic 
stresses.  This  analysis  will  include  deriving  the  equations 
governing  the  fluid  velocity  and  the  temperature  distribution. 

The  equations  governing  fluid  velocity  will  be  derived 
from  a  force  balance  approach.  We  shall  consider  the  forces 
that  act  on  a  differentially  small  parcel  of  fluid  to 
determine  its  behavior. 

The  equations  governing  temperature  will  be  derived  from 
the  principle  of  conservation  of  energy.  Energy  and 
temperature  are  closely  related.  In  fact,  in  an 
incompressible  fluid  temperature  is  a  direct  measurement  of 
internal  energy. 

These  equations  will  then  be  programmed  to  provide  a 
computer  simulation  for  predicted  velocity  and  temperature 
fields  for  various  parameters.  These  simulations  will  tell  us 
whether  or  not  it  is  possible  to  "shape"  velocity  and 
temperature  distributions  using  magnetic  fields.  Possible 
applications  include  heat  exchanges  and  any  transfer  process 
using  fluid  flow  as  a  transport  medium. 
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Abstract 


This  study  is  an  analysis  of  incompressible  flow  of  a  Newtonian 
fluid  past  a  vertical,  flat  plate  with  thermal  and  magnetic  fields. 
This  analysis  will  include  deriving  the  equations  governing  the 
fluid  velocity  and  the  temperature  distribution. 

The  equations  governing  fluid  velocity  will  be  derived  from  the 
conservation  of  linear  momentum.  Gravitational  forces,  thermal 
forces,  electromagnetic  forces,  and  viscous  forces  are  considered. 

The  equations  governing  temperature  will  be  derived  from  the 
principle  of  conservation  of  energy.  Energy  and  temperature  are 
closely  related.  In  fact,  in  an  incompressible  fluid  temperature 
is  a  direct  measurement  of  internal  energy. 

These  equations  will  then  be  programmed  to  provide  a  computer 
simulation  for  predicted  velocity  and  temperature  fields  for  various 
parameters.  These  simulations  will  tell  us  whether  or  not  it  is 
possible  to  "shape"  velocity  and  temperature  distributions  using 
magnetic  fields.  Possible  applications  include  heat  exchangers  and 
any  transfer  process  using  fluid  flow  as  a  transport  medium. 
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Chapter  One 


Definition  of 
Variables  and  the 
Geometry  of  the  Study 


The  first  step  in  formulating  a  solution  to  a  problem 
in  fluid  mechanics  is  to  understand  the  geometry  of  the 
problem  and  define  all  variables  involved. 

This  problem  involves  a  Newtonian  Fluid  in  a  steady, 
incompressible  flow  past  a  vertical  plate  with  thermal  and 
magnetic  fields. 


U_  ,  T„ 


Figure  1-1 

Note:  this  is  a  left-handed  system,  therefore  i  x  j  =  -k . 
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The  following  variables  will  be  used  throughout  the 
discussion : 

v  =  v  =  velocity  =  (u  I  +  v  j  +w  £) ,  — 

s 

u  =  velocity  in  the  x-direction,  — 

s 

v  =  velocity  in  the  y-direction,  — 

s 

w  =  velocity  in  the  z-dizection,  — 

s 

L  =  length  of  plate,  m 

T0  =  Temperature  of  plate ,  K 

U„  =  free  stream  velocity,  — 

s 

T„  =  free  stream  temperature,  K 

u  =  absolute  viscosity  =  0.001  -^2  for  water 

m-s 

p  =  density  =  1000  -^2  for  water 

mz 

2 

=  kinematic  viscosity  =  -^  =  l XI O'6  —  for  water 

P  s 


v 
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o  =  charge  density, 


rcr 


B  =  B  =  magnetic  field  strength, 


kg 

C'S 


k  =  coefficient  of  thermal  conductivity  , -J^L 

s*K 


C  =  heat  capacity  { constant  pressure ) , 


kg-K 


M  =  magnetic  index  number  = 


Re  =  Reynold's  number  = 

v 


LoB 

Ujp 


The  physical  situation  being  analyzed  is  an 
incompressible,  Newtonian  fluid  flowing  past  a  flat  plate  of 
length  L  at  a  rate  U„.  The  plate  is  at  uniform  temperature 
T0.  The  ambient  fluid  is  at  a  different  temperature  T_  . 
Define  the  x-direction  as  the  direction  along  the  length  of 
the  plate,  the  y-direction  as  the  direction  perpendicular  to 
the  plane  of  the  plate,  and  the  z-direction  as  the  vertical 
direction,  perpendicular  to  both  the  x-  and  y-direct ions . 
Additionally,  define  u  as  the  velocity  component  in  the  x- 
direction,  v  as  the  velocity  component  in  the  y-direction, 
and  w  as  the  velocit  component  in  the  z-direction. 

The  analysis  is  limited  to  two  dimensions  such  that  the 
velocity  field  is  a  function  of  x  and  y  and  the  temperature 
field  is  a  function  of  x  and  y. 
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Chapter  Two _ 

Derivation  of 
Equations  Governing 
Velocity  Distribution 

In  this  chapter,  the  partial  differential  equations 
governing  fluid  motion  are  derived  from  the  conservation  of 
linear  momentum.  The  result  is  two  equations,  one  for  the 
balance  of  forces  in  the  x-direction  and  one  for  the  balance 
of  forces  in  the  y-direction. 

These  equations  are  derived  from  Newton's  Second  Law  of 
motion : 


F=4  (jdv)  2-1 

at 

For  an  incompressible  flow,  this  means: 

E^PdV-f  2-2 

Where  dV  is  a  differential  element  of  volume. 

For  steady  two-dimensional  flow  Eq.  2-2  can  be  written  as: 


dy 


2-3 


Two  types  of  forces  are  considered:  body  forces  and 
surface  forces.  Body  forces  are  forces  that  may  be 
considered  to  act  through  the  center  of  mass  of  the  control 
volume  of  the  fluid.  The  body  forces  that  are  treated  are 
gravity,  thermal  forces,  and  electromagnetic  forces. 
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Surface  forces  act  on  the  surface  of  a  control  volume  of 
fluid.  Surface  forces  are  expressed  as  the  integral  of  the 
stresses  over  the  control  surface.1 

Gravity  exerts  a  force  on  a  different ially-small 
element  of  fluid  in  the  following  manner: 

Fg=pcNg  2-4 


Figure  2-1 

Because  this  force  acts  only  in  the  z-direction  one  neglects 
its  effects  in  horizontal  planar  flow. 

Thermal  forces  are  convective  and  are  caused  by  changes 
in  fluid  density  due  to  temperature2. 

Fc=PdVgp(r1-r2)  2-5 

Since  this  force  acts  in  the  z-direction  it  may  be 
neglected. 


Robert  A.  Granger,  Fluid  Mechanics  (New  York:  Holt, 
Rhinehart  and  Winston,  1985)  187. 

2Adrian  Bejan,  Convective  Heat  Transfer  (New  York:  John 
Wiley  and  Sons,  1984)  114. 


The  electromagnetic  force  is  expressed  as3: 


FB=gv  x  B=acNv  x  B 


2-6 


4  A  A  A  4  > 

JL 

B 


Figure  2-2 

Sigma  is  the  charge  density  and  B  is  the  magnetic  field 
strength.  Considering  only  forces  in  the  x-  and  y- 
directions  Eq.  2-6  simplifies  to: 

Fb=ocN(  -vBzi  +  uBzj)  2-7 

Consider  next  the  surface  forces.  Figure  2-3  shows  a 
dif ferent ially-small  element  of  fluid  with  dimensions 
(dx,dy,dz),  with  assorted  stresses  being  applied.4 


3Paul  A.  Tipler,  Physics  for  Scientists  and  Engineers  (New 
York:  Worth  Publishers  Inc.,  1991)  783. 

4Dr .  Hermann  Schlichting,  Boundary  Laver  Theory  (New  York: 
McGraw-Hill  Book  Co.  1968)  253. 
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°2Z+  5ldz 


Wnere : 

x^sheax  stxess  in  the  b-dixection  fiom  sheax  in  the  a-plane 

o  =  stxess,  NOT  chaxge  density 
Summing  the  forces  in  the  x-direction 

£  P dydz*  dz~x xz)  pdxdv 


+  ^xy+-5r Tdy-*xy)  Pdxdz 


y  i  a.  i ;  i  -i.  —  v—.. 


**  cy 
2-8  one  obtains 


2-8 


Tfs={  +^Z  +^£ )  pdxdvdz 
*  dx  dy  dz 

. t. i  1  a r  1 ,  the  forces  ir.  the  v— direction  simciifv  to 


2-9 


9 


£  V  <  *  "w 1  p  dxdyd* 


From  fluid  dynamics5 


/  dv  du \ 

Xxy=xyx~V,(-fa  +  -fy) 


x  =x  =u(i^+iT) 

y*  *y  *  dy  dz 


_  __  =11 ,  dw.  6uv 
xz  zx  3x 


2-10 


2-lla 


2-llb 


2-llc 


„  0U 


2-lid 


°r=-P*2P|^  2‘lle 

Substituting  Eq  2-lla  through  2-lle  into  Eq  2-9  yields: 


Y*  =[-^P+u(-^H+ZH)]dV 

^  s*  dx  V  5x2  3y2 


2-12 


Similarly,  the  surface  forces  in  the  y  direction,  Eq  2-10, 
simplify  to 


Yf =[-.§p+u(J^+i^)]dv 
^  sy  0y  3x2  dy2 


2-13 


Summing  the  force  components  in  the  x-direction  using  Eq  2- 
3,  2-7,  and  2-12  results  in: 


/ 

3x2 


+  .^H)  ]  dV-ovB  dV 
dy2 


2-14 


Since  the  plate  is  infinitesimally  thin,  the  pressure 


^Granger,  184. 


is  uniform  in  both  the  x-  and  y-directions,  hence  any 
derivatives  of  pressure  with  respect  to  x  and  y  are  zero. 
Using  this  information,  Eq  2-14  simplifies  to: 
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du  du  11  ,  d*u  d*u  \ 


2-15 


dx  dy  p  dx2  dy2  P 
The  equation  for  conservation  of  linear  momentum  in  the  x- 
direction  will  be  referred  to  as  Eq  A.  Similarly  the 
conservation  of  linear  momentum  in  the  y-direction  becomes: 


dv  dv  ll  ,  d^v  d^v  \  o£z 

dx  dy  p  dx2  dy2  P 


This  will  be  referred  to  as  Eq  B. 


2-16 
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Chapter  Three _ 

Formulation  of  the 
Energy  Equation 

The  First  Law  of  Thermodynamics  states  that  the  energy 
accumulated  in  a  control  volume  is  equal  to  the  energy 
entering  the  control  volume  m?nus  the  energy  leaving  the 
control  volume. 

The  energy  entering  the  control  volume  (C.V.)  and  the 


energy  leaving  the  control  volume  is  composed  of6: 

Rate  of  energy  accumulation  in  the  C.V.=  [1] 

Rate  of  transfer  of  energy  by  fluid  flow  +  [2] 

Rate  of  heat  transfer  by  conduction  +  [3] 

Rate  of  internal  heat  generation  -  [4] 

Rate  of  net  work  transfer  from  C.V.  to  the  [5] 

environment . 


These  terms  can  be  expressed  in  the  following  manner: 

[1]  =kxLy-~  (pe)  3-1 


Since  we  assume  the  flow  to  be  steady,  the  time  rate  of 
change  of  internal  energy  is  zero,  so  [1]  =  0. 

From  figure  3-1, 


[2]  =pveAx-  [pve+-^  (pve)  Ay]  Ax+pueAy-  [pue+-^  (pue)  Ax]  Ay  3-2 


6Bejan,  9. 


I  (  pve+  -j-  (  pve)dy)dx 


Figure  3-1 

q  3-2  simplifies  to 
[2]  =  -pA*A  y{u4^?+v|jf ) 

From  Figure  3-2  one  sees  that 


6qy 

+  T-* 


dy)dx 


[3]  =-  (AxAy)  +  =  -AxAyV-g 
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The  Fourier  Law  of  heat  transfer  is 
g=-JcVT  3-5 

Thus,  Eq.  3-4  becomes: 

[3]  =AxAy/cV2r  3-6 

From  Figure  3-2 

[4]  =AxAyg  3-7 

No  heat  is  generated  in  our  fluid,  so  q' ''=0  and  [4]  =  0. 


In  order  to  calculate  [5],  the  rate  of  net  work 
transfer  from  the  C.V.  to  the  environment,  the  forces  that 
exist  at  the  interface  between  the  control  volume  and  the 
environment  must  be  known.  As  defined  in  Chapter  Two,  these 
forces  are  called  surface  forces.  The  rate  of  work  transfer 
is : 

[5]=W=JLfFs-dr  3-8 

But  for  steady  flow,  the  surface  forces  do  not  change  with 
respect  to  time.  Thus,  Eq  3-8  simplifies  to: 

[5]=/fyi>  3-9 

The  above  term  is  analogous  to  heat  produced  by  fluid 
friction.  The  flows  that  are  of  concern  are  slow  enough 
that  turbulence  may  be  ignored.  Hence  the  work  done  by  the 
surface  forces  will  be  neglected.  In  order  to  determine  if 
this  is  a  realistic  assumption,  temperature  profiles 
predicted  by  the  resulting  energy  equation  should  be 
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compared  to  experimental  results. 

Combining  all  of  the  above  terms  results  in 

[1]  =  [2]  +  [3]  +  [4]  -  [5] 

0  =  -pCpAxAy(u~+v-^)  +AxAykV2T+0-Q 

Rewriting  the  above  yields 

„dT  dT__  k  ,  &T, 

dx  dy  p CpK  dx2  dy2 

This  equation  will  be  referred  to  as  Eq  C. 


3-10 


3-11 
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Chapter  Four _ 

Prandtl  Order  Reduction  and 

Conversion  to 
Dimensionless  Variables 

The  problem  has  been  reduced  to  three  equations  in 
terms  of  three  unknowns:  u,  v,  and  T  that  are  functions  of 
space,  x  and  y.  The  equations  may  be  simplified  somewhat  by 
determining  if  any  of  the  terms  are  relatively  insignificant 
and  may  be  neglected.  The  technique  of  Prandtl  Order 
Reduction  can  be  used  to  accomplish  this. 

The  first  step  in  the  reduction  process  is  to  assign 
each  variable  an  order  of  magnitude.  This  indicates  the 
size  of  the  values  that  the  particular  variable  might 
represent  in  comparison  to  other  variables.  For  example, 
the  velocity  u  may  be  assigned  an  order  of  magnitude  of  1 
and  the  velocity  v  an  order  of  magnitude  of  5,  because  the 
expected  values  of  v  are  much  smaller  than  the  expected 
values  of  u.  Similarly,  x  has  an  order  of  magnitude  of  1 
and  y  has  an  order  of  magnitude  of  8  because  previous  work 
done  in  fluid  mechanics  indicates  that  velocity  is 
essentially  constant  except  for  a  small  region  close  to  the 
plate,  so  one  expects  the  y-dimension  to  be  much  smaller 
than  the  x-dimension.  Recall  Eq  2-15 
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=-&  { .§!“  +  v 

dy  p  0x2  3y2  P 


Analyzing  the  orders  of  magnitude: 


2-15 


+6 


l=Ji(JL+A)-^ 

8  p  i2  62  p 


4-1 


On  the  left  side  of  the  equation,  all  terms  are  of  order  of 
magnitude  1.  The  two  viscous  terms,  however,  are  of 
different  orders  of  magnitude.  One  has  order  of  magnitude  1 
while  the  other  has  an  order  of  magnitude  1/82,  which  is 
much  larger  than  one.  Thus,  we  neglect  the  term  with  order 
of  magnitude  1.  Eq  2-15  becomes: 


u^H  +viH  =Ji  &E -22* v 
dx  dy  p  dy2  P 

Using  the  same  method,  Eq  2-16  becomes: 


4-2 


U!r*v4^=-*M— )♦ 

dx  dy  p  dy2 


d2v .  ^ 


u 


4-3 


There  now  exist  three  equations  relating  the  variables 
u,  v,  T,  x,  and  y.  To  write  these  equations  in  terms  of 
dimensionless  parameters,  define  the  following  dimensionless 
qualities : 
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4-4a 


4— 4b 


4-4c 


4-4d 


0= 


T-T0 

T~~T0 


4-4e 


Using  Eq  4-4a  through  4-4e,  Eq' s  4-2,  4-3,  and  3-11  are 
transformed  to  dimensionless  parameters. 

Recall  Eq  4-2 

4-2 

dx  dy  p  dy 2  p 


oB , 


Uj2  dx  dy  uj  P  dy 2  P 


v) 


4-5 


U  u.  +  V  U-  _ _ £ 


JJ 

a. 


*1 


3<l>2 


Las, 

P^« 


v 


4-6 


Eq  4-2  then  becomes: 


it -mv 

d£  dr|  Se  2 


4-7 


(Recall  our  definition  of  M  and  Re  from  Chapter  One.) 


Similarly,  Eq  4-3  becomes: 
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u-2£_*v-*£.-±?£-*uu- 

SC  Sri  Re  3t|2 

If  we  hold  To  constant 


4-8 


_3  ,r_T  x  _  dT  _  dT0  _  dT 
dx '  0  Sx  3x  Sx 

The  same  relationship  holds  for  the  derivative  of  (T-T0) 
with  respect  to  y  and  second  derivatives  with  respect  to 
both  x  and  y.  We  may  use  this  to  convert  Eq  3-11  to 
dimensionless  form. 


4-9 


(  — )  ( 

uj 


t-t 


X  ,  d(T-T0)  d(T~T0) 

)  (u - _2-+v - 2-)  = 


Sx 


3y 


Ls,  1  w_  k,,&(T-T0)  d2(T-TQ) 


(  — )  ( 
u 


T-T, 


)  (-  — )  ( 

pcD 


dx2 


dy 2 
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This  simplifies  to: 


T-T 
d-  0 


0L 


21- r. 


0.+.V 


r -r. 


o  ^  - 


si? 

X, 


CL 


si: 

L 


)  =- 


?CDUJ. 


{ 


S2  -  r°- 
T-Tn 


d{V 


s2-^ 

r.-r0 

3(^)2 

Lj 


4-11 


Thus,  Eq  3-11  becomes: 


u.se+v.se=__j__{Sf0  +  s^) 

SC  Stj  pCp£7^v5C2  Sri2 


4-12 


Equations  4-12,  4-7,  and  4-8  do  not  fully 
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mathematically  state  the  problem.  A  set  of  boundary 
conditions  must  be  included.  The  boundary  conditions  come 
from  physical  constraints  of  the  fluid  flow. 

The  first  constraint  is  called  the  "no  slip 
condition."7  This  states  that  the  velocity  of  the  fluid  at 
the  plate  is  zero. 

u*= 0  at  t|=0  BC  (!) 

v*=0  at  n=o 

The  second  constraint  is  that  at  distances  infinitely 
far  from  the  plate  the  velocity  approaches  the  free  stream 
velocity . 

v^U„i  as  t]  '°°  BC  (2) 

This  leads  to 

u*-l  as  T|-<®  BC  (2) 

The  third  boundary  condition  is  at  distances  infinitely 
far  from  the  plate  the  shear  stress  is  zero. 

(~  +-f^)  -o  as  y”°°  BC  (3) 

oy  ox 

The  technique  of  Prandtl  Order  Reduction  shows  that  the 
partial  derivative  of  v  with  respect  to  x  may  be  neglected: 

7Stuart  Churchill,  Viscous  Flows:  the  Practical  Use  of 
Theory  (Boston:  Butterworth  Publisher,  1988)  256. 
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t=m4+t)  4-13 

O  1 

The  boundary  condition  of  zero  shear  stress  simplifies  to 

—zr—  “*0  as  T)  "*0  BC  (3) 

dTj 

Similar  conditions  exist  for  the  temperature 
distribution,  T.  From  the  definition  of  the  function  0,  0 
must  equal  0  at  the  plate: 


0= 


T-T0  _  T0-T0 
T-Tn  T-Tr, 


=  0  at  tj=0 


BC  (4) 


Similarly  0  must  approach  1  at  distances  infinitely  far 
from  the  plate: 


0= 


r.-r0 

Tm~T0 


=i 


as  T|-o° 


BC  (5) 


A  final  boundary  condition  exists  on  the  temperature 
distribution.  At  distances  infinitely  far  from  the  plate, 
the  temperature  distribution  should  reach  equilibrium  and 
therefore  no  heat  should  be  transferred  This  means  that: 

»-0  as  tj-oo  BC  (6) 

drj 

These  boundary  conditions,  along  with  the  equations. 


constitute  a  sec  of  partial  differential  equations  which 
model  the  problem. 
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Chapter  Five 


Using  these  relationships.  Eg  4-7  becomes: 


ay  yr.aryr,  1  a3?*  ay 

dr\  dr)df  d(  a^2  Re  a^3  dC 

Sinilarly,  Eq  4-8  becomes: 


5-3 


5-4 


_  ay  a2?*  A  ay  a2?* .  1  a3?*  >  x„ay  5_5 

3t|  ac2  ac  OTiac  i?e  ari2ac  ^ 

Now  the  boundary  conditions  must  be  described  in  terms 
of  the  stream  function.  The  no  slip  condition  becomes: 
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9Y*  _  9T '  _a  i  «  Bp  /i  \ 

- ^-=0  at  t]  =0  BC  (1) 

The  boundary  condition  for  the  streamwise  velocity 
infinitely  far  from  the  plate  becomes: 

^-1  asti-»  BC  (2) 

dii 

The  boundary  condition  on  shear  stress  becomes: 

as  ri-oo  BC  (3) 

dx\2 

For  the  first  attempt  at  a  solution,  *F  was  assumed  to 
have  the  form: 

(B£crid) 

The  partial  derivatives  of  that  appear  in  Eq  5-4  are: 


U^AbCV^U)  +ABdCa*Vd-1f/<€) 

dt] 


5-6a 


-^^=Ab{i3-l)CV2^(€)  +ABd(2Jb+d-l)Ca*cTiw‘1f/(€)  + 

dti2 

AB2d2^a*2c^b*2d-2fH^  5-6b 

^l=Ab{b- 1)  (b-2)  Cari b~3f(e)  + 

atl3 

ABd[3b(b+2d-2)  +{d- 2)  (d-1)  Ca+ctlw'3/'(e)  + 

AB23d2  (jb+d-2)  £a+2CllJ>-2d-3.f//(€)  +AB3d3Ca*3cTii’*3d‘3/ (3)  (e)  5-6c 

-Sjjp  =Aaf a"  (e)  +  ABcf  a+c-1ii  b*df'  { c )  5-6d 

^2m* 

3.-,  -AabQa~xx\b~x  f  (e)  +AS(jbc+ad+ad) 

d^dri 

+AB2cdCa*2c-1rii>+2d-1f//(€)  5-6e 

Where,  for  notational  purposes,  define  e  such  that: 

€=BCctlef  5-7 

Substituting  Eq  5-6a  through  5-6e  into  Eq  5-4  and 


simplifying 
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AabZa~xr\b~2  f2  (e)  + 

AB(bc+ad+bcd-ad2)  (e)  /'(e)  + 

AB2d(c+ad-bc)  ft {e)  f/(c)  + 

AB2d(bc+ad)  {•*2c-iT,i»*2<f-2jf(C)  _f//{€) 

_i3(i?-l)  (i?~2)  ^r1^_3f(€)  -AMaC^r,^ (€) 

_  Bd [3-b (i?+2d~2  )  +  (d~l)  (d-2)  ]  ^a*c^i+d-3  f  /  -ABMcCa+c~1Tl't>+c,// (e) 

Rg 

_  3B2cf(d-l)  (d~2)  ^a*2Cj^b*2d-3 f  //  _  B2d2  ^a-*-3c^2)*3d-3^- (3)  _q5_Q 

i?e 

The  similarity  transformation  was  chosen  to  transform 
the  partial  differential  equation  involving  the  function  *F 
into  an  ordinary  differential  equation  involving  the 
function  f.  An  ordinary  differential  equation  is  one  where 
all  of  the  coefficients  of  the  derivatives  contain  only 
constants  and  powers  of  the  argument  of  the  function  f. 

This  means  that  a,  b,  c,  and  d  must  be  chosen  in  such  a  way 
that  the  coefficients  of  the  derivatives  of  f  only  involve 
powers  of  e,  but  do  NOT  contain  any  free  powers  of  £  or  f\. 
In  order  to  search  for  values  of  a,  b,  c,  and  d  that  will 
satisfy  these  requirements,  a  table  was  created  that  lists 
the  coefficient  of  each  term  and  the  exponents  of  £  and  r\ 
for  that  term,  as  well  as  the  exponents  after  powers  of  e 


were  factored  out. 


Using  this  definition  of  4**,  u’  and  v*  become 
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u*=Ai\  [e-f'U)  +2jf(e)  ]  5-10 

vn=AB£~2r\sfl(e)  5-11 

The  resulting  boundary  conditions  are: 

u*=Ar\  [e/'U)  +2f  (c)  ]  =0  at  tj=0  (1) 

v*=AB('2ti5.f'(€)  =0  at  Tj  =0  BC  (*) 

The  no  slip  conditions  are  met. 

u*=A<x>[<x>f/ (<*>)  +2f(<*>)  ]  =1  at  r\=<*>  (2) 


This  cannot  be  used  as  a  boundary  condition  for  the  computer 
solution.  The  T]  which  appears  in  Eq  5-10  makes  it 
impossible  for  u*  to  satisfy  the  boundary  conditions.  If 
one  defines  f(°°)=f'  (<»)  =0,  one  obtains  u"  =  0  as  i)  approaches 


infinity,  not  u*=l  as  T\  approaches  infinity,  which  is  a 
physical  constraint  of  any  solution  found. 
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Chapter  Six _ 

Second  Solution  Attempt: 

Blasius  Analysis 

The  first  similarity  transformation  attempted  did  not 
give  a  viable  solution.  The  next  step  is  to  examine  a 
solution  to  a  similar  problem  that  had  been  solved  many 
years  ago.  The  solution  to  the  problem  of  flow  past  a  flat 
plate  without  thermal  or  magnetic  fields  was  solved  in  1902 
by  Blasius  using  the  method  of  selecting  a  similarity 
transformation  to  change  the  partial  differential  equation 
of  the  conservation  of  linear  momentum  in  the  x-direction 
into  an  ordinary  differential  equation.8 

Blasius  assumed  the  following: 

U=C7.F{€)  ;  v=[7„G(e)  6-1 

where 


and  5  was  a  function  of  x  alone. 

The  streamline  function  becomes 

y  e 

JP=fudy=U'„j’F{e)  50cfe  =  U.60f  (e)  6-3 

0  0 


where 


8Granger,  713. 


6-4 


f{€)  -Jf(€)  dc6-4 


This  definition  yields 


v=-||*-  lu.S'of  <*>  *U.S0f'(e)  -*] 


6-6 


v*  =  -6n/(e)  +6q  ef'(e) 


6-7 


6-8 


3u*_  02^‘  _  f  "(e) 


di\  dr\2  5C 


6-9 


au:=afY:=_^€jf//(€) 

ac  acarj  60er 


6-10 


a2u,.63y,.fw(€) 
at)2  aq3 


6-11 


Define : 


*'-d6o.  f/fr>_cff(€) 
60  —  ,  f  (e)  gg- 

Recall  Eq  4-7 


*au*  *  du*  i  a2u*  „  * 

U  +  V  =  -±-  -  Wv 

3C  ari  i?e  0tj2 

Substituting  6-6  through  6-11  into  Eq  4-7 
Sofi'of  (e)  /"(€)  (€)  +Mt>20b'0efl(e)  =0 

KG 


6-12 


4-7 


6-13 


In  order  to  make  this  an  ordinary  differential  equation,  all 
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of  the  coefficients  of  the  derivatives  must  involve  only 
constants  and  powers  of  e.  For  this  to  be  true  two 
conditions  must  be  satisfied  at  once: 


60 


d50 

dC 


/ 


6-14 


The  only  solution  that  satisfies  this  condition  is  for  80  to 
be  a  constant.  This  means  that  epsilon  is  a  function  of  y 
alone,  which  in  turn  means  that  u,  which  is  a  function  of 
epsilon,  depends  on  y  alone.  Experience  dictates,  however, 
that  u  depends  also  on  x. 
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Chapter  Seven 


Third  Solution  Attempt: 
Modified  Blasius  Method 

The  third  solution  attempt  involves  slightly  modifying 
the  Blasius  method.  Instead  of  assuming: 


u=tr„F(e)  ;  e=-£ 

o„ 


u  was  assumed  to  have  a  slightly  different  form 


u=U„xaybF(e)  ;  €  =  -jr—  7 

5o 

Following  an  analysis  similar  to  that  done  in  the  Blasius 
solution  for  a  flat  plate  presented  in  Chapter  Six 


J  J 

xF=Judy=j U„xaybF(-^~)  dy=UegxaJ'  {eb0)bF{e)  6 0de 


x¥=U'jcatf*1f(z)  ;  f  (e)  =  f  ebF(e)  de 


Using  Eq  7-4,  one  obtains  the  following  expressions: 


u=-~=U„xa5$f,(e) 


7-5a 


v=~~~  =  ~Uj<a~1^o*1f  (e)  -U.(i)+l)xa5Xf(€)  +£J^a5Xef'(e) 


7 -5b 


=  ^^  =  C/“aXa"15°f/(€)  +C7-faxa5o'lfiof/<€)  7-5c 
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JiU|* 


7-5d 


&u_  a3^ 
dy2  dy3 


=C7-xa5f'1f///(€) 


7-5e 


Entering  Eq  7-5a  through  7-5e  into  Eq  4-2  yields  the 
following: 


Ulax^-'bff'te)  /'(e)  +  cdbx2aa2i'1a£fr/(e)  f'(e) 


-Ofx^a^aoe/^e)  /"(e)  -ulax2*-xbf£{e)  f"{e) 

+  Ufx2a5f_16/0€f/(e)  /"(e)  =  2.Ujxab$'2£"(e)  +j®°  [/..ax^a^1/ (e) 

P  P 


+^0.(^1)  xaafaU(e)  -^iuaaX€f'(e)  7-6 

p  p 

Simplifying  Eq  7-6  by  dividing  by  U„  xa  80b-2  yields: 
UAax^b^bx^*1^]  f'f'-Ujax a-1a§+2+  (^lJx^X]  iff" 

J£ax-xa^f+^£  (i?+i)  a^of-^a^a^f^o  7-7 

p  p  p  p 

A  80  function  that  makes  this  an  ordinary  differential 
equation  must  now  be  found.  There  are  a  number  of 
coefficients  that  involve  8's  and  powers  of  x: 

1)  xa*18o *2  =  Ci;  2)  xa8o+18o  =  C2;  3)  x^a^CJ;  4)  a2ao  =  C4 
Working  with  condition  (3)  gives: 

i 

x^a^CJ  -  a0  =  C3x1 


7-8 


This  is  interesting  because  it  indicates  that: 


which  is  similar  to  the  transformation  that  changed  the 
equation  from  a  partial  differential  equation  into  an 
ordinary  differential  equation  in  Chapter  Five.  In  that 
transformation,  a  function  of  y3  over  x  was  used,  whereas  in 
this  solution  a  function  of  y  over  x1/3  is  used. 

Next  one  must  check  to  see  if  this  function  of  80 
satisfies  the  other  conditions.  Looking  at  condition  4: 

_i  .2 

5q6o=(C3*  3)2(-|C3x  3)  =±03*= Constant  1~10 

Thus  condition  4  is  satisfied. 

Looking  at  condition  1 

J.  3a-3  *b*2 

xa'18o+2=xa'1  (C3x  3  )b*2  =  Cx  3  =Constant  -  3a+£>-l  =  0 

Looking  at  condition  2 

_1  __2  3a*^>*l-2 

xa6o*18o=xa(CJx  3)i>+1(-|cJx  3 )  =Cx  3  - 

=  Constant  -  3a+b-l  =  0  7-12 

Thus  3a  +  b  -  1  =  0  is  a  relationship  between  a  and  b 
that  must  exist  in  order  to  satisfy  conditions  1  and  2.  Now 
examine  the  boundary  conditions.  Recall  from  Eq  7-2  that 
this  analysis  began  with  the  assumption: 


7-11 
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u=UmxaybF(e)  -  u*=xaybF(e)  7-2 

One  of  the  boundary  conditions  is  that  u  approaches  1 
as  y  approaches  infinity  (BC  2) .  But  as  y  approaches 
infinity,  u  approaches  infinity  for  this  definition  of  u. 

One  way  to  attempt  to  solve  this  dilemma  is  to  modify  the 
boundary  condition  to  state  that  u*  approaches  1  as  y 
approaches  some  multiple  of  80  rather  than  allowing  y  to 
approach  infinity. 

_1  3  3  a+b 

u*=xaybF{-¥T)=xa(Cx3)bF(±Z-r)=cbx  3  F(C)  =1  -  3a+b=0  7-13 

x1  x 1 

But  3a+b=0  and  3a+b=l  are  inconsistent  equations.  A 
solution  that  matches  the  boundary  conditions  does  not 
exist . 

Note  that  this  is  the  same  situation  encountered  with 
the  first  solution  in  Chapter  Five.  In  the  first  solution  a 
transformation  was  found  that  yielded  an  ordinary 
differential  equation,  but  could  not  satisfy  the  boundary 
conditions.  Here,  using  a  very  similar  transformation,  the 
same  problem  is  encountered. 
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Chapter  Eight _ 

Fourth  Solution  Attempt: 
Constrained  Magnetic  Field 

Recall  from  the  second  solution  attempt  in  Chapter  Six 
utilizing  the  Blasius  method  of  searching  for  a  similarity 
transformation  that  Eg  6-13  was: 


4-f (3)  (€)  +Mo f  (€)  f"U)  -Mb2Xf{e)  +M&2Xef'{e)  =0 

KG 

Where  we  had  defined  the  following  in  Eq  6-12 


5o  =  ^;  f'(€) 
d(. 

Let  M  =  M0TTl. 


_  dfU) 
de 

Recall : 


6-13 


6-12 


c-JL  -  60--3- 

60  € 


8-1 


Using  Eq  8-1  and  the  new  definition  of  M  to  rewrite  Eq  6-13: 

-Lfiv  (e)  +MU(e)  f"U)  ~M0 rr^Xf  (€)  +A^tf1-3 -606/0ef/(e)  =0  8-2 
Re  €  € 


Simplifying  Eq  8-2  yields 


4-f (3)  (e)  +Mo^(e)  f"(*)  -W0-60 b'Qf(e)  +Af0506^f/(c)  =0  8-3 

Re  6 

Now  one  must  change  this  into  an  ordinary  differential 
equation.  Examining  the  coefficients  of  the  derivatives  in 
Eq  8-3,  one  finds  that  50  80'  must  be  a  constant. 

Mo=c  8-4 

80  =  0  at  x  =  0 .  Thus,  Eq  8-4  yields 
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60-v^S;  S'® 

Verifying  Eq  8-5  satisfies  Eq  8-4 

Mo=(W  (1^)  =  (1)  (2C)  -C 

Recall  that  for  the  transformation  the  following 
expressions  were  formulated  (See  chapter  Six) : 

u=U„£'ie)  6-8 

v=-uXfU)+uXef'(e)=(Uj  (±x~?)  [ef'(e) -f  (e)  ]  6-6 

Using  the  boundary  conditions  on  u  and  v  to  establish 
the  boundary  conditions  on  the  function  f (E) . 

u=U„f'(e)  =0  at  y=0  -  f'(0)  =0  BC  (1) 

v»(-|)  (Ujc  2 )  [ef'U)  -Jf(€)  ]  =0  at  y=0  -  f(0)  =0  BC  (!) 

u=U„f/(e)  =  U„  as  y-~,  =1  BC  (2) 

This  constitutes  an  ordinary  differential  equation  with 
a  set  of  boundary  conditions  that  may  be  solved  using  a 
computer.  Using  an  approach  known  as  the  "shooting"  method, 
Professor  Malek-Madani  of  the  United  States  Naval  Academy 
Math  Department  wrote  a  program  for  Mathematical  software 
that  was  capable  of  solving  these  equations.  The  solution 
method  involves  treating  the  problem  as  an  initial  value 
problem  and  iteratively  converging  the  solution  to  the 


36 


underlying  boundary  value  problem.  See  appendix  A  for  a 
copy  of  the  Mathematical  program. 
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Chapter  Nine _ 

Discussion  of  Results 

The  computer  simulation  was  used  to  see  what  variations 
in  the  magnetic  field  did  to  the  velocity  profile.  The 
simulation  was  run  and  the  outputs  are  shown  for  various 
"Magnetic  Index"  numbers  and  Reynolds  Numbers.  These 
numbers  are  a  representation  of  the  strength  of  the  magnetic 
force  and  the  viscous  force  exerted  upon  the  fluid  flow. 

Recall  from  Chapter  One  that  the  definition  of  the 
Magnetic  Index  number  was: 


Some  typical  figures  are 
L  =  25  cm 

a  =  0.001  Molarity  (.001  mole  free  charges  /  Liter) 

Bz  =  1.0  Gauss  (Approximately  twice  the  strength  of  the 
Earth's  magnetic  Field.) 

U  =  7  cm/s 
p  =  1000  kg/m3. 

These  figures  yield  a  Magnetic  Index  number  M  =  0.0344  and  a 
Reynolds  Number  of  Re  =  17,500  .  This  gives  an  idea  of  the 
magnitude  of  the  numbers  that  should  be  used  in  the  computer 
simulation.  Appendix  B  contains  computer  simulations  for 
various  Re  and  M  values . 

These  graphs  include  plots  of  u*  and  v*  versus  position 
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for  various  Reynolds  numbers  and  Magnetic  Index  numbers. 
Three  values  of  Reynolds  number  were  used:  13000,  17500,  and 
22000;  for  each  of  these  Reynolds  numbers  three  values  of 
the  Magnetic  Index  number  were  used:  0  (to  show  the  velocity 
profile  with  no  magnetic  field),  0.5,  and  -0.5  .  One  trial 
with  M  =  1.5  was  included  (Figure  B-ll)  to  show  exaggerated 
effects  of  the  magnetic  field.  For  the  simulation  with  Re  = 
17500  and  Mo  =  0.5  graphs  of  f (£) ,  f' (£)  and  f ' ' (£)  versus 
epsilon  were  also  included. 

These  graphs  are  believable  because  they  show  a  number 
of  characteristics  the  equations  developed  from  fluid  flow 
theory  indicate  should  appear. 

First  of  all  it  can  be  seen  that  the  "boundary  layer", 
the  line  where  the  streamwise  velocity  approaches  the  free 
stream  velocity,  closely  follows  the  boundary  layer 
predicted  by  the  classical  Blasius  Solution  for  flow  without 
a  magnetic  field.  This  can  be  seen  by  comparing  any  of  the 
simulation  outputs  with  a  graph  of  the  boundary  layer 
predicted  by  classical  theory, 

°  predicted 

A  graph  of  the  predicted  boundary  layer  versus  £  is  included 
after  the  graph  of  the  first  simulation  output  as  Figure  B- 
3.  Compare  Figure  B-l  and  Figure  B-3  and  you  can  see  that 
for  zero  magnetic  field,  this  simulation  predicts  a  boundary 
layer  at  T\  equals  approximately  0.035  for  £  =  1.  Figure 
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B-3,  the  predicted  boundary  layer  according  to  the  Blasius 
Solution,  shows  a  boundary  layer  of  0.037  for  £  =  1. 

A  second  detail  shown  in  the  figures  that  is  supported 
by  theory  is  that  the  velocity  perpendicular  to  the  plate, 
v,  is  affected  by  changes  in  M  to  a  larger  degree  than  the 
streamwise  velocity,  u.  Compare  Figure  B-l  and  B-2  with 
Figure  B-4  and  B-5  to  see  this  difference.  For  the  same 
change  in  M,  the  u  distribution  does  not  change  much  (Figure 
B-l  compared  to  Figure  B-4) .  A  slight  hump  is  produced  at 
the  boundary  layer  in  Figure  B-4.  The  v  distribution, 
however,  changes  much  more  noticeably  (Figure  B-2  compared 
to  Figure  B-5) .  In  Figure  B-5  we  see  that  the  v 
distribution  peaks  at  a  value  of  0.015  and  settles  to  a 
value  of  approximately  0.005,  as  compared  to  Figure  B-2, 
where  the  v  distribution  peaks  at  a  value  of  0.02  and 
settles  there.  This  difference  in  the  sensitivity  to 
changes  in  the  magnetic  field  makes  sense,  because  the 
magnetic  force  involves  the  cross  product  of  the  magnetic 
field  and  the  velocity  of  the  charged  particles  (See  Eq  2- 
6) .  The  velocity  in  the  streamwise  direction  (x-direction) 
is  much  larger  than  the  velocity  perpendicular  to  the  plate 
(y-direction) .  This  means  that  the  magnetic  force  operating 
in  the  y-direction  is  greater  than  the  magnetic  force 
operating  in  the  x-direction,  so  a  change  in  M  should  have  a 
greater  affect  on  the  velocity  in  the  y-direction. 

Along  the  same  lines,  we  notice  that  the  "hump" 
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produced  by  M  =  0.5  is  slightly  larger  for  larger  Reynolds 
Numbers.  This  can  be  seen  by  looking  at  Figure  B-15  and 
Figure  B-21.  Figure  B-15  corresponds  to  a  slower  flow  than 
Figure  B-21,  and  the  "hump"  caused  by  the  magnetic  field  is 
less  noticeable  for  this  slower  flow.  This  corresponds  to 
theory  because  Reynolds  Number  is  directly  proportional  to 
the  free  stream  velocity.  A  higher  Reynolds  number  equates 
to  a  higher  velocity,  which  in  turn  equates  to  larger 
magnetic  forces. 

These  figures  show  some  interesting  trends.  One 
observation  that  can  be  made  is  that  the  slope  of  the 
velocity  fields,  both  u  and  v,  is  affected  by  changing  M. 

By  comparing  Figures  B-l  and  B-4,  one  can  see  that  as  M 
increases  from  0  to  0.5,  u  approaches  the  free  stream 
velocity  more  quickly.  In  fact,  u  appears  to  overshoot  the 
free  stream  velocity  and  then  settle  back  to  match  the  free 
stream  velocity.  This  overshoot  is  even  more  exaggerated  in 
Figure  B-ll,  where  u  reaches  a  value  of  1.2  times  the  free 
stream  velocity  then  settles  to  a  final  value  matching  the 
free  stream  velocity.  Conversely,  u  appears  to  "flatten 
out"  when  M  changes  from  0  to  -0.5,  as  can  be  seen  by 
comparing  Figures  B-l  and  B-9.  The  velocity  u  reaches  a 
higher  final  value  in  Figure  B-9  than  in  Figure  B-l,  but 
approaches  this  final  velocity  more  slowly. 

The  effect  the  of  magnetic  field  on  the  slope  is  even 
more  evident  in  the  graphs  of  the  v  distributions. 
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Comparing  Figures  B-2  and  B-5  show  that  as  the  magnetic 
field  increases  from  0  to  0.5,  v  changes  the  same  way  that  u 
changed.  Figure  B-5  shows  that  the  v  distribution  rises 
more  sharply  than  in  figure  B-2,  then  settles.  Similarly, 
Figure  B-10  shows  that  for  a  negative  M,  the  v  distribution 
rises  more  slowly  but  settles  at  a  higher  value  than  the  v 
distribution  in  Figure  B-2.  This  indicates  that  magnetic 
fields  can  be  used  to  distribute  flows  more  evenly  or  create 
sharper  boundary  layers. 

The  figures  also  show  that  negative  values  of  M  result 
in  higher  velocities  at  greater  distances  from  the  plate. 
This  can  be  seen  by  comparing  figure  B-l  and  B-9,  B-13  and 
B-17,  and  B-l 9  and  B-23.  In  each  case  the  velocities  u  and 
v  reach  higher  values  for  negative  values  of  M.  It  may  be 
possible  to  use  magnetic  fields  to  control  the  speed  of 
flows . 

What  these  simulations  indicate  is  that  the  velocity 
profile  can  be  shaped  by  using  a  magnetic  field.  This  could 
have  many  applications.  Fluid  flow  is  used  mostly  as  a 
process  of  transporting  something,  whether  it  is  heat  being 
transported  by  a  pipe  in  a  radiator  or  oxygen  being 
transported  by  blood  in  an  artery,  and  these  processes  of 
transportation  often  depend  on  fluid  velocity.  Thus, 
controlling  the  shape  of  the  velocity  profile,  allows 
control  of  the  transport  process.  One  could  transport  heat 
through  a  pipe  with  less  heat  loss  if  one  could  shape  the 
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velocity  profile  properly.  One  could  aid  in  the  diffusion 
of  oxygen  across  artery  membranes  by  controlling  the 
boundary  layer  of  blood  flow. 

The  next  step  that  needs  to  be  taken  in  this  research 
is  to  validate  the  computer  simulations  by  running  trials 
and  taking  measurements  to  verify  the  predicted  results  with 
experimental  data.  This  study  was  supposed  to  generate 
empirical  data,  but  failed  to  do  so.  This  was  due  in  part 
to  the  fact  that  the  laser  doppler  velocimeter  was 
inoperable.  A  hot  film  anemometer  was  tried,  but  this 
caused  a  number  of  problems. 

In  addition  to  this,  the  computer  simulation  must  be 
expanded  to  include  Eq  (C) ,  the  Energy  Conservation 
Equation,  so  a  graph  of  the  predicted  temperature  profile 
can  be  created  in  the  same  way  that  graphs  of  u  and  v  were 
created. 
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Htxl _ x2_,  x3_,  tj  -  x2;  g[x1_,  x2_,  x3_,  t_]  -  x3; 

re  =  17500;  MO  =  .5; 

Cl  =  re  MO/2; 

C2  *  re/2; 

h[x1_,  x2_,  x3_,  tj  =  Cl  t*(-1)  xl  -  Cl  x2  -  C2  xl  x3; 

(“•We  define  f«f  (epsilon),  g*f  '(epsilon),  and  h  ■  f"  (epsilon)*) 

(■"Then  we  solve  our  equation  for  f'" (epsilon)  =  h'.*) 

tstart  ■  0.000001; 
tfinal  -  6; 

F[b_]  :=  Module  [{sol}  ,  sol  =  NDSol  ve[{x1 '  [t]  ==  f[x1[t],  x2  [t] ,  x3[t],t], 

x2'[t]  ==  g [xl  [t3 ,  x2[t],  x3[t]  ,t] , 
x3'[t]  ==  h [xl  [t] ,  x2 [t] ,  x3 [t] ,t] , 

xl [tstart]  *»  0,  x2[tstart]  ==  0,  x3[tstart]  ==  b},  {xl ,  x2,  x3}, 
{t,  tstart,  tfinal}]; 

out  =  Evaluate  [{xl  [t] ,  x2[t],  x3[t]}  /.  sol/,  t  ->  tfinal]; 

{out  [[1,2]]-  1}] 

(This  module  takes  has  the  initial  value  of  f" (epsilon)  as  a  variable*) 

(*and  gives  the  final  value  of  f' (epsilon)  -  1  as  output*) 

e  =  Fi ndRoot [F [b]  [ [1 ] ] ,  {b,  1,  0.S}]; 

(This  tells  us  which  initial  value  of  f"  minimizes  *) 

(*the  value  of  f' (tfinal)  -  1  *) 

soil  *  NDSol ve [{xl ' [t]  =a  f  [xl  [t],  x2[t],  x3[t],t], 
x2'  [t]  ««  g[x1[t],  x2[tl ,  x3[t]  ,t] , 
x3'[t]  -*  h[x1  [t] ,  x2[t] ,  x3[t],t], 

xl [tstart]  «*  0,  x2[tstart]  ««  0,  x3[tstart]  *=  e[[l,2]]}, 

{xl,  x2,  x3} ,  {t,  tstart,  tfinal}]; 

(This  solves  the  system  of  equations  we  have  set  up  using*) 

(*the  initial  value  for  f"  we  calculated  using  the  FindRoot  function.*) 

plotl  *  PlotCxI  [t]  /.  soil,  {t,  tstart,  tfinal}, 

AxesLabel  ->  {“epsilon", "f (epsilon)"}, 

PlotLabel->{"  -  Re"  re,"  -  Mo,  graph  of  f(z)“M0}] 

PSPrint  [plotl] 

plot2  =  Plot[x2[t]  /.  soil,  {t,  tstart,  tfinal), 

AxesLabel  ->  {"epsilon", “f '(epsilon)"} , 

PlotLabel->{"  *  Re"  re,"  *  Mo,  graph  of  f'(z)"M0}3 

PSPrint[plot2] 

pi ct3  =  Plot[x3[t]  /.  soil,  {t,  tstart,  tfinal}, 

AxesLabel  ->  {"epsilon", "f "(epsilon)"} , 

Plot  Label ->{“  -  Re"  re,"  =  Mo,  graph  of  f"(z)"M0}] 


PSPrint [pi ot3] 


piot3  =  Plot [x3[t]  /.  soil,  {t,  tstart,  tfinal}, 

AxesLabel  ->  {"epsilon", "f"(epsilon)"}, 

PlotLabel->{"  *  Re"  re,"  *  Mo,  graph  of  f"(z)"M0}] 

PSPrint[plot3] 

uttj  =  x2[t]  /.  soil; 

v[tj  =  .5  (t  x2[t]  -  xl  [t])  /.  soil; 

plot4  -  Pi ot3D  [u [y/Sqrt [x] ]  [ [1 3  3 ,  {x,  0.1,  1}.  {y,  0.00001,  .4}, 

PlotLabel  ->  {"  «  Re"  re,  "  *  Mo*  MO, "graph  of  u(x,y)"}, 
AxesLabel  ->  rzeta","etaV'u(x,y)"}, 

PlotPoints->  30] 

pi otB  -  Plot3D[xA(-.5)  v[xA(-.5)  y][[l33,  {x,  0.1,  1},  {y,  0.00001,  .4}, 
PlotLabel  ->  ‘graph  of  vCx.y)". 

AxesLabel  ->  {"zeta","eta","v(x,y)'}, 

PlotPoints->  30] 


Simulation  Trials 


Figure  B-9 


graph  of  v(x,y) 
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Figure  B-10 
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graph  of  v<x,y> 


Figure  B-14 


Figure  B-18 


